
REAL(8) FUNCTION func_gN(x0)
    
USE GLOBAL
USE OMP_LIB

implicit none

REAL(8),INTENT(IN)::x0(6)

REAL(8)::v2(na,nb2),c1(na,nb2),uc1(na,nb2)
REAL(8)::wmax,p0max,h0max,gn0max,ct0max,m0max,tau0max,ww0max,f0max,b_next_global
REAL(8)::value1,value2,value3
INTEGER  i_b, i_a, i_default_global,now_time
INTEGER::index_b,new_index_b,number_of_threads,i_b_base,try
INTEGER(4)::hours2,mins2,zero_h
INTEGER(4)::IRIGHT,ILEFT
REAL(8)::a_initial,b_initial,udef,u_value,indicator_external,gg
REAL(8)::stats(6),xlb(6),xub(6),penalty,time,ct_matrix_new(nb, na)
EXTERNAL udef
INTEGER(4), PARAMETER::nb3=501
REAL(8)::end_time_loop,start_time_loop      
REAL(8)::objective_u_fe,objective_u_unemp
EXTERNAL objective_u_fe,objective_u_unemp
REAL(8)::objective_u_fe_gN,objective_u_unemp_gN,objective_fun_gN_wbar,objective_u_fe_wbargN
EXTERNAL objective_u_fe_gN,objective_u_unemp_gN,objective_fun_gN_wbar,objective_u_fe_wbargN
REAL(8)::bgrid_base(nb3),h_baseX(na,nb3),gn_baseX(na,nb3),ct_baseX(na,nb3),bp_baseX(na,nb3),pn_baseX(na,nb3)    
REAL(8)::w_base(na,nb),pn_base(na,nb),ct_base(na,nb),aux0
REAL(8),DIMENSION(nb,na)::v_matrix_new,v0_matrix_new,vaut_matrix_new
REAL(8)::p0,term1,cnval,ctval 
REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb),DCSVAL,epsilongrid(ncdf), dev_q_nodef
INTEGER::i0(1),i
REAL(8)::value_rule,value_optimize
REAL(8)::x_matrix_new2(nb,na),gn_matrix_new2(nb,na),ct_matrix_new2(nb,na)

do i=1,nb2
   bgrid2(i) = bmin + (bmax - bmin) * REAL((i-1d0) / (nb2 - 1d0))
end do

i0 = LOCATE2(bgrid2,0d0)
i_bzero2 = i0(1)
bgrid2(i_bzero2)=0d0


index_bfeas2(:) = 0
iterfun = iterfun+1

call compute_bgrid

beta =x0(1)
chi  =x0(2)
icost=x0(3)
psi1 =x0(4)
psi2 =x0(5)
wbar =x0(6)

!Initialize terminal node with default as risk-free debt environment
q_matrix_nodef = payoff/ (lambda + r)
v0_matrix   = 0d0
vaut_matrix = 0d0

do i_a = 1,na         
    FDATA(1:nb) = v0_matrix(:,i_a)
	ILEFT  = 0
	IRIGHT = 0
	DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(2) - bgrid(1))
	DRIGHT = (FDATA(nb)-FDATA(nb-1)) / (bgrid(nb) - bgrid(nb-1))
	call  DCSDEC (nb, bgrid, FDATA(1:nb), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, 1:nb), coeff_matrix(i_a, :,1:nb))
ENDDO

do i_a = 1,na
    FDATA = q_matrix_nodef(:,i_a)
	ILEFT  = 0
	IRIGHT = 0
    DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(2) - bgrid(1))
	DRIGHT = (FDATA(nb)-FDATA(nb-1)) / (bgrid(nb) - bgrid(nb-1))
	call  DCSDEC (nb, bgrid, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, BREAK_GRID, CSCOEF)
    break_matrix_q(i_a, :) = BREAK_GRID
    coeff_matrix_q(i_a, :,:) = CSCOEF
end do

    
indicator_external =1  !FROM EXTERNAL FILE
if (indicator_external < 0.5) then
    
        open (10, FILE='graphs\v.txt',STATUS='replace')
        open (11, FILE='graphs\default.txt',STATUS='replace')
        open (12, FILE='graphs\x.txt',STATUS='replace')
        open (13, FILE='graphs\gn.txt',STATUS='replace')
        open (14, FILE='graphs\b_next.txt',STATUS='replace')
        open (16, FILE='graphs\ctt.txt',STATUS='replace')

        do i_b = 1,nb 
            do i_a = 1,na 

                a_initial = agrid(i_a)
                b_initial = bgrid(i_b)
                try=1
                i_default_global = 2 !defaults  

                IF (choice_constant_gn.EQ.1)  constant_gn = constant_gn_value
                
                IF (psi1.LE.500d0) THEN
                    value1 = objective_u_fe(0d0,a_initial,b_initial,i_default_global)
                    value2 = objective_u_unemp(0d0,a_initial,b_initial,i_default_global)                                
                
                    u_value = MAX(value1,value2)             
                    vaut_matrix_new(i_b, i_a) = u_value-udef(a_initial)                    
                 ELSE
                    vaut_matrix_new(i_b, i_a) = -1d10                    
                 ENDIF
                 
                i_default_global = 1 !repays                 
                b_next_global=(1d0-lambda)*b_initial
                
                IF (choice_constant_gn.EQ.1)  constant_gn = constant_gn_value+coefyt*a_initial
                
                b_next_global=bmin
               
                value1 = objective_u_fe_gN(b_next_global,constant_gn,a_initial,b_initial,i_default_global)
                value2 = objective_u_unemp_gN(b_next_global,constant_gn,a_initial,b_initial,i_default_global)
                value3 = objective_u_fe_wbargN(b_next_global,constant_gn,a_initial,b_initial,i_default_global)
                value_rule = MAX(value1,MAX(value2,value3))

                IF (value1.GE.MAX(value2,value3)) THEN
                    x_matrix_new(i_b, i_a)=x_global1
                    gn_matrix_new(i_b, i_a)=gn_global1
                    ct_matrix_new(i_b, i_a)=ct_global1
                ELSEIF (value2.GE.MAX(value1,value3)) THEN
                    x_matrix_new(i_b, i_a)=x_global2
                    gn_matrix_new(i_b, i_a)=gn_global2
                    ct_matrix_new(i_b, i_a)=ct_global2
                ELSE
                    x_matrix_new(i_b, i_a) =x_global3
                    gn_matrix_new(i_b, i_a)=gn_global3
                    ct_matrix_new(i_b, i_a)=ct_global3
                ENDIF               

                value_optimize = -1d11
                
               gn_matrix_new2(i_b, i_a)=1d3
               
                IF (gn_matrix_new2(i_b, i_a).LE.constant_gn_value) THEN                    
                   u_value = value_optimize
                   x_matrix_new(i_b, i_a)=x_matrix_new2(i_b, i_a)
                   gn_matrix_new(i_b, i_a)=gn_matrix_new2(i_b, i_a)
                   ct_matrix_new(i_b, i_a)=ct_matrix_new2(i_b, i_a)
                ELSE
                    u_value = value_rule
                ENDIF

                v0_matrix_new(i_b, i_a) = u_value  
                v_matrix_new(i_b, i_a) = MAX(vaut_matrix_new(i_b, i_a), v0_matrix_new(i_b,i_a))

                 if (vaut_matrix_new(i_b, i_a) <= v0_matrix_new(i_b, i_a)) then
                   default_decision(i_b, i_a) = 1
                    b_next_matrix(i_b, i_a) = b_next_global
                 else
                   default_decision(i_b, i_a) = 2
                    b_next_matrix(i_b, i_a) = 0d0
                 end if              

          
                 WRITE(10, '(F18.10, X, F18.10, X, F15.10)') MAX(-100d0,v_matrix_new(i_b, i_a)),MAX(-100d0,v0_matrix_new(i_b, i_a)),MAX(-100d0,vaut_matrix_new(i_b,i_a))
                 WRITE(11, '(F18.10)') defaultgrid(default_decision(i_b, i_a))
                 WRITE(12, '(F18.10)') x_matrix_new(i_b, i_a)
                 WRITE(13, '(F18.10)') gn_matrix_new(i_b, i_a)
                 WRITE(14, '(F15.11)') b_next_matrix(i_b, i_a)
                 WRITE(16, '(F18.10)') ct_matrix_new(i_b, i_a)
                 
                 IF (psi1.LT.500d0) THEN        !default vs. risk-free debt economy         
                     q_matrix(i_b, i_a) = 0d0                 
                     q_matrix_nodef(i_b, i_a) = 0d0 
                 ELSE
                     q_matrix(i_b, i_a) = payoff/ (lambda + r)
                     q_matrix_nodef(i_b, i_a) = payoff/ (lambda + r)
                 ENDIF
                 
            end do
      end do
  
    CLOSE(10)
    CLOSE(11)
    CLOSE(12)
    CLOSE(13)
    CLOSE(14)
    CLOSE(16)

    v_matrix    = v_matrix_new
    v0_matrix   = v0_matrix_new
    vaut_matrix = vaut_matrix_new

                
    PRINT*, 'No Initial Guess Uploaded'
    PRINT*,' '
else !READ DATA FROM EXTERNAL FILES

open (10, FILE='graphs\v.txt')
open (16, FILE='graphs\q_paid.txt')
open (17, FILE='graphs\b_next.txt')
open (23, FILE='graphs\bounds.txt')

      READ(23, '(F15.11, X, F15.11)') bmin, bmax

      do i_b = 1,nb
         do i_a = 1,na

             READ(10, '(F18.10, X, F18.10, X, F18.10)') v_matrix(i_b, i_a), &
                     v0_matrix(i_b, i_a), vaut_matrix(i_b, i_a)
             READ(16, '(F15.11, X, F15.11)') q_matrix(i_b, i_a), q_matrix_nodef(i_b, i_a)
             READ(17, '(F15.11)') b_next_matrix(i_b, i_a)
             if (vaut_matrix(i_b, i_a) > v0_matrix(i_b, i_a)) then
                 default_decision(i_b, i_a) =2
             else
                 default_decision(i_b, i_a) =1
             end if 

         end do
      end do
CLOSE(10)
CLOSE(16)
CLOSE(17)
CLOSE(23)

PRINT*, 'Saved Initial Guess Uploaded'

end if

!       beta      chi     icost     psi1    psi2     wbar
xlb = [0.900d0,  0.002d0,  -0.1d0,   0.1d0,   0.5d0,     2.4d0]
xub = [0.945d0,  0.85d0,  0.5d0,   0.7d0,   9.5d0,   3.4d0]

call iterate_gN

call pol_fcns_gN 

!CALL INITIAL_gN

call MC_SIM(stats)

func_gN = ((-stats(1)-22.8d0)/22.8d0)**2d0 +((stats(2)-18.1d0)/18.1d0)**2d0 +((stats(3)-1.1d0)/1.1d0)**2d0 +((stats(4)-1.4d0)/1.4d0)**2d0+((stats(5)-1.11d0)/1.11d0)**2d0+((stats(6)-0.9815d0)/0.9815d0)**2d0

call system_clock(count=now_time)      
time  = (now_time-start_time)/rate_time/MAX(1,iterfun-1)
hours2 = INT(FLOOR(time/3600))
mins2  = FLOOR((time/3600-(hours2))*60d0)
PRINT*,'Iter = ',iterfun, 'Time per iteration',hours2,'hs',mins2,'mins'

PRINT*, 'MODEL SOLVED AND MC STATS COMPUTED'
RETURN

END FUNCTION func_gN
!******************************************************************************************
REAL(8) FUNCTION objective_u_fe_gN(x,gg,a_initial2,b_initial2,i_default_global2)

USE GLOBAL

IMPLICIT NONE

REAL(8),INTENT(IN)::x,gg,a_initial2,b_initial2
INTEGER,INTENT(IN)::i_default_global2
REAL(8)::pres,hres,gnres,cnres,ctres,cres,a0,b0,aux0
REAL(8)::wres,taures,ares,wwres,utilc,utilg,gnval,hval
REAL(8)::bpval,wval,tauval,ct_global
REAL(8)::p0,m0,h0,q,b
REAL(8)::ctval,cnval,cval,q_fun,penalty,fval
REAL(8)::aux_lhs,balt,aux0alt,minc,maxc,maxc2,maxc3
REAL(8)::term1,taxbar,bound1,bound2,new_value,old_value
REAL(8)::objective_u_feold,gn_global1old,x_global1old,ct_global1old 
REAL(8)::q2,aux02,p02,gnval2 ,temp1,temp2
INTEGER,PARAMETER::num=200
REAL(8)::ctemp(num)
INTEGER(4)::i,jj,im,iax,feas,error,niter2,ftype,try,index_opt,i_num
LOGICAL::flag00
EXTERNAL q_fun

objective_u_feold = -1d10

b0 = b_initial2
a0 = DEXP(a_initial2)

b = x

q = q_fun(b,a_initial2)

aux0 = -q*(b-(1d0-lambda)*b0) + payoff*b0           !p*gn - T
aux0 = aux0*(1d0-defaultgrid(i_default_global2))

ctval= a0+aux0
ct_global = ctval   
IF (ctval.LE.0d0) THEN
    objective_u_fe_gN= huge_n +1d3*(0d0-ctval)             
    if (choice_glob.EQ.1) THEN
        gn_global1 = 0d0
        x_global1  = 0d0
        ct_global1 = 0d0
    ENDIF
    RETURN
ENDIF

gnval = gg
cnval = 1d0-gnval
ctval = a0+aux0

p0   = (1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)
wval = alpha*p0
tauval = p0*gnval-aux0

taxbar = taxrate*(a0+in*p0)

flag00 = ((p0.LE.0d0).OR.(wval.LT.wbar).OR.(gnval.GT.1d0).OR.(ctval.LE.0d0).OR.(tauval.GT.taxbar)) 

IF (.NOT. flag00) THEN
    pres  = p0
    hres  = 1d0
    ctres = ctval
    taures= tauval
    wwres = wval
    gnres = gnval
    feas  = 1d0
ELSE
    ctres   = 0d0    
    feas    = 0d0
    taures  = tauval
    hres    = 0d0
    pres    = p0
    gnres   = gnval
    wwres   = wval
ENDIF

!negative penalty
penalty = 1d3*((pres.LE.0d0)*(0d0-pres)+(wwres.LT.wbar)*(wbar-wwres)+(ctres.LE.0d0)*(0d0-ctres))
IF (taxrate.LE.1d0) penalty = penalty +1d6*((taures.GT.taxbar)*(taures-taxbar))
penalty = penalty+1d9*((b.LT.bmin)*(bmin-b))+1d2*((b.LT.bmin))
penalty = penalty +1d3*((gnres.GT.1d0)*(gnres-1d0)+(gnres.LT.0d0)*(0d0-gnres))

cnres = MAX(1d-8,hres**alpha-gnres)

IF (mu .NE. 0d0) THEN
    cres  = (omegac*ctres**(-mu)+ (1d0-omegac)*cnres**(-mu))**(-1d0/mu)
ELSE
    cres  = (ctres**(omegac))*(cnres**(1d0-omegac))
ENDIF

cres = cres - (cnres.LE.1d-6)*(-1d6)
cres = cres - (ctres.LE.1d-6)*(-1d6)

IF (sigma.EQ.1d0) THEN
    utilc = DLOG(MAX(cres,1d-8))
ELSE 
    utilc = MAX(cres,1d-8)**(1d0-sigma)/(1d0-sigma) 
ENDIF

IF (sigmag.EQ.1d0) THEN 
    utilg = DLOG(MAX(gnres,1d-8))
ELSE
    utilg = MAX(gnres,1d-8)**(1d0-sigmag)/(1d0-sigmag) 
ENDIF

wres = (1d0-chi)*utilc+chi*utilg 
wres = feas*wres + (1d0-feas)*huge_n +penalty 

if (choice_glob.EQ.1) THEN
    gn_global1 = gnres
    x_global1  = 1d0
    ct_global1 = ctres
ENDIF

objective_u_fe_gN = wres

end FUNCTION objective_u_fe_gN
!***************************************************************************
REAL(8) FUNCTION objective_u_unemp_gN(x,gg,a_initial2,b_initial2,i_default_global2)

USE GLOBAL

IMPLICIT NONE

REAL(8),INTENT(IN)::x,gg,a_initial2,b_initial2
INTEGER,INTENT(IN)::i_default_global2
REAL(8)::pres,hres,gnres,cnres,ctres,cres,a0,b0,aux0
REAL(8)::wres,taures,ares,wwres,utilc,utilg,gnval,hval
REAL(8)::lambda1,pp1,bpval,wval,tauval
REAL(8)::p0,m0,h0,q,b,x_global2old,objective_u_unempold 
INTEGER,PARAMETER::num=200
REAL(8)::htemp(num),aahat,aa1    
REAL(8)::ctval,cnval,cval,q_fun,penalty,fval,gn_global2old,ct_global,ct_global2old
REAL(8)::aux_lhs,balt,qalt,aux0alt,minh,maxh,temp1,taxbar,netbc,bound1,bound2,new_value,old_value
INTEGER(4)::i,jj,iax,feas ,error,niter2,try,ftype,i_num,index_opt
REAL(8)::objective_u_unemp_gNold
INTEGER(4)::i_num_star,i_num_star2,index_optmax,root
real(8)::sign_new,sign_old,sign_two,ctemp2,difc,old_valuemax
LOGICAL::flag00
EXTERNAL q_fun

objective_u_unempold = -1d10
root=0

b0 = b_initial2
a0 = DEXP(a_initial2)
b = x

q  = q_fun(b,a_initial2)

aux0 = -q*(b-(1d0-lambda)*b0) + payoff*b0            !p*gn - T
aux0 = aux0*(1d0-defaultgrid(i_default_global2))
aux_lhs   = aux0+taxrate*(a0+in*wbar/alpha)                          !always smaller than cT
!aux_lhs < 0 follows from  pn*gn = aux0+tax < 0 
!                          alpha*p*h**alpha=wbar*h < wbar

IF (choice_nowbar.EQ.1) aux_lhs = 100d0     !no aux_lhs condition for wbar=0

ctval  = a0+aux0
ct_global = ctval

IF( (choice_fe.EQ.1).OR. (choice_nowbar.EQ.1)) THEN
    objective_u_unemp_gN= huge_n
    if (choice_glob.EQ.1) THEN
        gn_global2= 0d0    
        x_global2 = 0d0
        ct_global2= 0d0
    ENDIF
    RETURN
ELSEIF  (aux_lhs.LE.0d0) THEN
    balt = b+1d-6
    q    = q_fun(balt,a_initial2) 
    aux0alt = -q*(balt-(1d0-lambda)*b0) + payoff*b0            
    aux0alt = aux0alt*(1d0-defaultgrid(i_default_global2))
    objective_u_unemp_gN= huge_n+1d3*(aux0alt-aux0)             
    if (choice_glob.EQ.1) THEN
        gn_global2 = 0d0    
        x_global2  = 0d0
        ct_global2 = 0d0
    ENDIF
    RETURN    
ENDIF

!case 1: h<1 and nonbinding tax cap
try = 1
gnval = gg

ftype = 199
minh  = (gnval+1d-4)**(1d0/alpha)
!Alternative: minc = (alpha/wbar*(1d0-omegac)/omegac)**(1d0/(1d0+mu))*ctval           !b8
maxh  = 1d0

old_value = AFunction(ftype,a0,gnval,ct_global,minh)
old_valuemax = 10d+33
i_num_star=0

index_opt = 0
index_optmax = 0
difc = (maxh-minh)/(REAL(num)-1d0)
       
140 DO i_num=1+i_num_star,num
    ftype = 199
    htemp(i_num) = minh+(maxh-minh)*(REAL(i_num)-1d0)/(REAL(num)-1d0)     
    new_value = AFunction(ftype,a0,gnval,ct_global,htemp(i_num))
    sign_new = SIGN(1d0,new_value)
    sign_old = SIGN(1d0,old_value)
    sign_two = SIGN(1d0,sign_new*sign_old)

    if (DABS(new_value) <= DABS(old_valuemax)) then
        index_optmax = i_num
        old_valuemax = new_value        
    end if      
        
    i_num_star = i_num
    if (sign_two.EQ.-1d0) then
        index_opt = i_num
        old_value = new_value
        GOTO 282
    ENDIF
    old_value = new_value

ENDDO       

282 IF (index_opt.EQ.0) THEN
    bound1=htemp(MAX(index_optmax-1,1))
    bound2=htemp(MIN(index_optmax+1,num))
ELSE
    bound1=htemp(MAX(index_opt-1,1))
    bound2=htemp(MIN(index_opt,num))        
ENDIF

hval=BrentRoots(ftype,a0,gnval,ct_global,bound1,bound2,1d-10,1000,fval,niter2,error)  

IF ((error.EQ.1)) THEN 
    hval = 0d0
    IF (i_num_star.LT.num)  go to 140
ELSEIF (error.EQ.2) THEN
    PRINT*,'Problem: More Iterations!'
    PAUSE
    STOP
ENDIF

cnval = hval**alpha-gnval
cnval = MAX(1d-6,cnval)
p0    = (1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)

1053 tauval = p0*gnval-aux0   
taxbar = taxrate*(a0+in*p0*hval**alpha)
    
flag00 = ((gnval.LT.0d0).OR.(gnval.GE.hval**alpha).OR.(p0.LE.0d0).OR.(ctval.LE.0d0).OR.(hval.GT.1d0).OR.(tauval.GT.taxbar)) 

IF (.NOT.flag00) THEN
    pres  = p0
    hres  = hval
    gnres = gnval
    ctres = ctval
    taures= tauval
    wwres = wbar
    feas  = 1d0
ELSE
    ctres = 0d0 
    feas  = 0d0
    hres  = 0.0d0
    gnres = 0d0
    taures= tauval
    pres  = p0
    wwres = wbar    
ENDIF
  
!negative penalty
penalty = 1d3*((pres.LE.0d0)*(0d0-pres)+(ctres.LE.0d0)*(0d0-ctres))
IF (taxrate.LE.1d0) penalty = penalty+1d6*((taures.GT.taxbar)*(taures-taxbar))
penalty = penalty+1d9*((b.LT.bmin)*(bmin-b))+1d2*((b.LT.bmin))
penalty = penalty +1d3*((gnres.GT.hres**alpha)*(gnres-hres**alpha)+(gnres.LT.0d0)*(0d0-gnres))
penalty = penalty +1d3*((hres.GT.1d0)*(hres-1d0)+(hres.LT.0d0)*(0d0-hres))

cnres = MAX(1d-8,hres**alpha-gnres)

IF (mu .NE. 0d0) THEN
    cres = (omegac*ctres**(-mu)+ (1d0-omegac)*cnres**(-mu))**(-1d0/mu)
ELSE
    cres = (ctres**(omegac))*(cnres**(1d0-omegac))
ENDIF

IF ((cnres.LE.1d-6).OR.(ctres.LE.1d-6)) cres = 0d0

aahat  = (hres+(1d0-hres)*alphac)**(-1d0)
aa1    = hres+(1d0-hres)*(alphac**(1d0-sigma))

IF (sigma.EQ.1d0) THEN
    utilc = (aahat**(1d0-sigma))*aa1*DLOG(MAX(cres,1d-8))
ELSE    
    utilc = (aahat**(1d0-sigma))*aa1*MAX(cres,1d-8)**(1d0-sigma)/(1d0-sigma)
ENDIF

IF (sigmag.EQ.1d0) THEN
    utilg = DLOG(MAX(gnres,1d-8))
ELSE
    utilg = MAX(gnres,1d-8)**(1d0-sigmag)/(1d0-sigmag)
ENDIF

wres = (1d0-chi)*utilc+chi*utilg 
wres = feas*wres + (1d0-feas)*huge_n + penalty 

objective_u_unemp_gN = wres
if (choice_glob.EQ.1) THEN
    gn_global2 = gnres
    x_global2 = hres
    ct_global2 = ctres
ENDIF

IF (root.EQ.0) THEN
    x_global2old  = x_global2
    gn_global2old = gn_global2
    ct_global2old = ct_global2
    
    objective_u_unemp_gNold = objective_u_unemp_gN
    root=1
ELSEIF (root.GT.0) THEN
    IF (objective_u_unemp_gN.GT.objective_u_unemp_gNold) THEN
        objective_u_unemp_gNold= objective_u_unemp_gN
        
        if (choice_glob.EQ.1) THEN
            gn_global2old = gn_global2    
            x_global2old  = x_global2 
            ct_global2old = ct_global2
        ENDIF
    ENDIF
    root=root+1
ENDIF
    
IF (i_num_star.LT.num)     GO TO 140

IF (objective_u_unemp_gNold.GT.objective_u_unemp_gN) THEN
    objective_u_unemp_gN= objective_u_unemp_gNold
        
    if (choice_glob.EQ.1) THEN
        gn_global2 = gn_global2old    
        x_global2  = x_global2old 
        ct_global2 = ct_global2old
    ENDIF
ENDIF
    

RETURN

end function objective_u_unemp_gN
!***************************************************************************
REAL(8) FUNCTION objective_u_fe_wbargN(x,gg,a_initial2,b_initial2,i_default_global2)

USE GLOBAL

IMPLICIT NONE

REAL(8),INTENT(IN)::x,gg,a_initial2,b_initial2
INTEGER,INTENT(IN)::i_default_global2
REAL(8)::pres,hres,gnres,cnres,ctres,cres,a0,b0,aux0
REAL(8)::wres,taures,ares,wwres,utilc,utilg,gnval,hval
REAL(8)::lambda1,pp1,bpval,wval,tauval
REAL(8)::p0,m0,h0,q,b,x_global2old,objective_u_unempold 
INTEGER,PARAMETER::num=200
REAL(8)::htemp(num),aahat,aa1    
REAL(8)::ctval,cnval,cval,q_fun,penalty,fval,gn_global2old,ct_global,ct_global2old
REAL(8)::aux_lhs,balt,qalt,aux0alt,minh,maxh,temp1,taxbar,netbc,bound1,bound2,new_value,old_value
INTEGER(4)::i,jj,iax,feas ,error,niter2,try,ftype,i_num,index_opt
LOGICAL::flag00
EXTERNAL q_fun

objective_u_unempold = -1d10

b0 = b_initial2
a0 = DEXP(a_initial2)
b = x

q = q_fun(b,a_initial2)

aux0 = -q*(b-(1d0-lambda)*b0) + payoff*b0            !p*gn - T
aux0 = aux0*(1d0-defaultgrid(i_default_global2))
aux_lhs   = aux0+taxrate*(a0+in*wbar/alpha)                          

IF (choice_nowbar.EQ.1) aux_lhs = 100d0     !no aux_lhs condition for wbar=0

ctval  = a0+aux0
ct_global = ctval

IF (choice_nowbar.EQ.1) THEN
    objective_u_fe_wbargN= huge_n
    if (choice_glob.EQ.1) THEN
        gn_global2= 0d0    
        x_global2 = 0d0
        ct_global2= 0d0
    ENDIF
    RETURN
ELSEIF  (aux_lhs.LE.0d0) THEN
    balt = b+1d-6
    q = q_fun(balt,a_initial2)

    aux0alt = -q*(balt-(1d0-lambda)*b0) + payoff*b0            
    aux0alt = aux0alt*(1d0-defaultgrid(i_default_global2))
    objective_u_fe_wbargN= huge_n+1d3*(aux0alt-aux0)          
    if (choice_glob.EQ.1) THEN
        gn_global2 = 0d0    
        x_global2  = 0d0
        ct_global2 = 0d0
    ENDIF
    RETURN    
ENDIF

gnval = gg
hval  = 1d0
p0    = wbar/alpha
cnval = 1d0-gnval
    
tauval = p0*gnval-aux0   
taxbar = taxrate*(a0+in*p0*hval**alpha)

IF (DABS(p0-(1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)).GT.1d-6) p0=0d0

flag00 = ((gnval.LT.0d0).OR.(gnval.GE.hval**alpha).OR.(p0.LE.0d0).OR.(ctval.LE.0d0).OR.(hval.GT.1d0).OR.(tauval.GT.taxbar)) 

IF (.NOT.flag00) THEN
    pres  = p0
    hres  = hval
    gnres = gnval
    ctres = ctval
    taures= tauval
    wwres = wbar
    feas  = 1d0
ELSE
    ctres = 0d0 
    feas  = 0d0
    hres  = 0.0d0
    gnres = 0d0
    taures= tauval
    pres  = p0
    wwres = wbar    
ENDIF
  
!negative penalty
penalty = 1d3*((pres.LE.0d0)*(0d0-pres)+(ctres.LE.0d0)*(0d0-ctres))
IF (taxrate.LE.1d0) penalty = penalty+1d6*((taures.GT.taxbar)*(taures-taxbar))
penalty = penalty+1d9*((b.LT.bmin)*(bmin-b))+1d2*((b.LT.bmin))
penalty = penalty +1d3*((gnres.GT.1d0)*(gnres-1d0)+(gnres.LT.0d0)*(0d0-gnres))
penalty=penalty+ 1d3*(wbar/alpha.GT.1d-6+(1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu))*(wbar/alpha - 1d-6-(1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu))
penalty=penalty+ 1d3*(wbar/alpha+1d-6.LT.(1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu))*((1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)-wbar/alpha-1d-6)

cnres = MAX(1d-8,hres**alpha-gnres)

IF (mu .NE. 0d0) THEN
    cres = (omegac*ctres**(-mu)+ (1d0-omegac)*cnres**(-mu))**(-1d0/mu)
ELSE
    cres = (ctres**(omegac))*(cnres**(1d0-omegac))
ENDIF

IF ((cnres.LE.1d-6).OR.(ctres.LE.1d-6)) cres = 0d0

aahat  = (hres+(1d0-hres)*alphac)**(-1d0)
aa1    = hres+(1d0-hres)*(alphac**(1d0-sigma))

IF (sigma.EQ.1d0) THEN
    utilc = (aahat**(1d0-sigma))*aa1*DLOG(MAX(cres,1d-8))
ELSE    
    utilc = (aahat**(1d0-sigma))*aa1*MAX(cres,1d-8)**(1d0-sigma)/(1d0-sigma)
ENDIF

IF (sigmag.EQ.1d0) THEN
    utilg = DLOG(MAX(gnres,1d-8))
ELSE
    utilg = MAX(gnres,1d-8)**(1d0-sigmag)/(1d0-sigmag)
ENDIF

wres = (1d0-chi)*utilc+chi*utilg 
wres = feas*wres + (1d0-feas)*huge_n + penalty 

objective_u_fe_wbargN = wres
if (choice_glob.EQ.1) THEN
    gn_global3 = gnres
    x_global3 = hres
    ct_global3 = ctres
ENDIF  

end function objective_u_fe_wbargN
!******************************************************************************************
!Compute the objective function in the Belman equation
DOUBLE PRECISION function objective_fun_gN(b,a_initial2,b_initial2,i_default_global2)

USE GLOBAL

IMPLICIT NONE

DOUBLE PRECISION,INTENT(IN)::b,a_initial2,b_initial2
INTEGER,INTENT(IN)::i_default_global2
INTEGER :: other_type, i,j,h, t, d, NINTV
REAL(8)::acum,exp_v_next,value_next,a_next,a_fun,scalar,DCSVAL,bb
REAL(8)::a_threshold,a_max1,a_min1,a_next1,acum2,acum1,g,q,DNORDF,v0_fun,vaut_fun
REAL(8)::gg,acum_int,acum_int_lo,acum_int_hi,borrowing, interpolate,cdf,cdf1,cdf_threshold, DNORIN,objective_u_unemp_gN
REAL(8)::exp_v_next_aut,acum3,objective_u_fe_gN
EXTERNAL DNORDF,a_fun,interpolate,v0_fun,vaut_fun,DCSVAL,DNORIN,udef,objective_u_unemp_gN,objective_u_fe_gN

REAL(8)::value1,value2 ,u_value,ucost,udef

IF (psi1.LT.500d0) THEN
    a_threshold   = a_fun(b,a_initial2)
ELSE
    a_threshold   = -1d4
ENDIF

cdf_threshold = DNORDF((a_threshold - rhoa*a_initial2 - (1-rhoa)* mua)/sigmaa)
bb = b*(1d0-defaultgrid(i_default_global2))

exp_v_next = 0d0
exp_v_next_aut = 0d0
scalar = 1d0 / SQRT(pi_number)

a_min1 = (1-rhoa)*mua + rhoa* a_initial2 - 4*sigmaa
a_max1 = (1-rhoa)*mua + rhoa* a_initial2 + 4*sigmaa

NINTV = ncdf - 1
acum1=0d+0
acum2=0d+0

if (a_threshold < a_min1) then

    do i=1,nquad_v
        cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
        a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
        value_next = v0_fun(bb,a_next)
        acum1 = acum1 + 0.5*(cdfmax-cdfmin)*quad_w_v(i) * value_next
    end do

elseif(a_threshold > a_max1) then

    do i=1,neps
        cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
        a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
        value_next = vaut_fun(a_next)
        acum1 = acum1 + 0.5*(cdfmax - cdfmin)*quad_w_v(i) * value_next
    end do

else
    do i=1,neps

    cdf = 0.5*(quad_x_v(i) + 1)*(cdf_threshold - cdfmin) + cdfmin
    a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
    value_next = vaut_fun(a_next)
    acum1 = acum1 + 0.5*(cdf_threshold - cdfmin)*quad_w_v(i) * value_next

!        compute integral for a > a_threshold
    cdf1 = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdf_threshold) + cdf_threshold
    a_next1 = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf1) * sigmaa 
    value_next = v0_fun(bb,a_next1)
    acum2 = acum2 + 0.5*(cdfmax - cdf_threshold)*quad_w_v(i) * value_next

    end do
end if

exp_v_next = (acum1+acum2)/(cdfmax - cdfmin)

IF (i_default_global2.EQ.2) THEN
    acum3=0d+0
        do i=1,neps
            cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
            a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
            value_next = vaut_fun(a_next)
            acum3 = acum3 + 0.5*(cdfmax - cdfmin)*quad_w_v(i) * value_next
        end do
    
        exp_v_next_aut = acum3/(cdfmax - cdfmin)
        exp_v_next = phi*exp_v_next + (1d0-phi)*exp_v_next_aut
ENDIF

ucost = defaultgrid(i_default_global2)*udef(a_initial2)

value1 = objective_u_fe_gN(bb,constant_gn,a_initial2,b_initial2,i_default_global2)
value2 = objective_u_unemp_gN(bb,constant_gn,a_initial2,b_initial2,i_default_global2)

u_value = MAX(value1,value2)


objective_fun_gN = -(u_value - ucost) - beta * exp_v_next
  
if (choice_glob.EQ.1) THEN
    IF (value1.GE.value2) THEN
        reg_global = 1
        x_global    = x_global1
        gn_global    = gn_global1
        ctrad_global = ct_global1  
    ELSE
        reg_global = 2
        x_global     = x_global2
        gn_global    = gn_global2
        ctrad_global = ct_global2  
    ENDIF    
ENDIF

end FUNCTION objective_fun_gN
!******************************************************************************************
!Compute the objective function in the Belman equation
DOUBLE PRECISION function objective_fun_gN_wbar(b,a_initial2,b_initial2,i_default_global2)

USE GLOBAL

IMPLICIT NONE

DOUBLE PRECISION,INTENT(IN)::b,a_initial2,b_initial2
INTEGER,INTENT(IN)::i_default_global2
INTEGER :: other_type, i,j,h, t, d, NINTV

REAL(8)::acum,exp_v_next,value_next,a_next,a_fun,scalar,DCSVAL,bb
REAL(8)::a_threshold,a_max1,a_min1,a_next1,acum2,acum1,g,q,DNORDF,v0_fun,vaut_fun
REAL(8)::gg,acum_int,acum_int_lo,acum_int_hi,borrowing, interpolate,cdf,cdf1,cdf_threshold, DNORIN,objective_u_fe_wbargN
REAL(8)::exp_v_next_aut,acum3
EXTERNAL DNORDF,a_fun,interpolate,v0_fun,vaut_fun,DCSVAL,DNORIN,udef,objective_u_fe_wbargN

REAL(8)::value1,value2 ,u_value,ucost,udef

IF (psi1.LT.500d0) THEN
    a_threshold   = a_fun(b,a_initial2)
ELSE
    a_threshold   = -1d4
ENDIF

cdf_threshold = DNORDF((a_threshold - rhoa*a_initial2 - (1-rhoa)* mua)/sigmaa)
bb = b*(1d0-defaultgrid(i_default_global2))

exp_v_next = 0d0
exp_v_next_aut = 0d0
scalar = 1d0 / SQRT(pi_number)

a_min1 = (1-rhoa)*mua + rhoa* a_initial2 - 4*sigmaa
a_max1 = (1-rhoa)*mua + rhoa* a_initial2 + 4*sigmaa

NINTV = ncdf - 1
acum1=0d+0
acum2=0d+0

if (a_threshold < a_min1) then

    do i=1,nquad_v
        cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
        a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
        value_next = v0_fun(bb,a_next)
        acum1 = acum1 + 0.5*(cdfmax-cdfmin)*quad_w_v(i) * value_next
    end do

elseif(a_threshold > a_max1) then

    do i=1,neps
        cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
        a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
        value_next = vaut_fun(a_next)
        acum1 = acum1 + 0.5*(cdfmax - cdfmin)*quad_w_v(i) * value_next
    end do

else
    do i=1,neps

    cdf = 0.5*(quad_x_v(i) + 1)*(cdf_threshold - cdfmin) + cdfmin
    a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
    value_next = vaut_fun(a_next)
    acum1 = acum1 + 0.5*(cdf_threshold - cdfmin)*quad_w_v(i) * value_next

!        compute integral for a > a_threshold
    cdf1 = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdf_threshold) + cdf_threshold
    a_next1 = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf1) * sigmaa 
    value_next = v0_fun(bb,a_next1)
    acum2 = acum2 + 0.5*(cdfmax - cdf_threshold)*quad_w_v(i) * value_next

    end do
end if

exp_v_next = (acum1+acum2)/(cdfmax - cdfmin)

IF (i_default_global2.EQ.2) THEN
    acum3=0d+0
        do i=1,neps
            cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
            a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
            value_next = vaut_fun(a_next)
            acum3 = acum3 + 0.5*(cdfmax - cdfmin)*quad_w_v(i) * value_next
        end do
    
        exp_v_next_aut = acum3/(cdfmax - cdfmin)
        exp_v_next = phi*exp_v_next + (1d0-phi)*exp_v_next_aut
ENDIF

ucost = defaultgrid(i_default_global2)*udef(a_initial2)

u_value = objective_u_fe_wbargN(bb,constant_gn,a_initial2,b_initial2,i_default_global2)

objective_fun_gN_wbar = -(u_value - ucost) - beta * exp_v_next

end FUNCTION objective_fun_gN_wbar
!********************************************************************************************
subroutine optimize_gN(b_next,v_value,a_initial2,b_initial2,i_default_global2)

USE GLOBAL

implicit none

REAL(8),INTENT(IN)::a_initial2,b_initial2
INTEGER,INTENT(IN)::i_default_global2
integer :: MAXFN, i,t, d, index_opt,index_opt2, index_vector(1), num
parameter (num=300)
DOUBLE PRECISION :: b_next, v_value, objective_fun_gN, new_value,  old_value,old_value2, b, vector(nb),&
                    b_next_grid(num), STEP, BOUND, XACC, b_next_guess
real(8)::b_next_inf,b_next_sup 
EXTERNAL objective_fun_gN, DUVMIF
REAL(8)::b_next2,v_value2,bmax_aux

REAL(8)::objective_u_unemp_gN,objective_u_fe_gN,value3,value1,value2
EXTERNAL::objective_u_unemp_gN,objective_u_fe_gN

IF (taxrate.LT.1d0) THEN
    b_next_inf = b_initial2-0.2d0
    b_next_sup = b_initial2+0.2d0
ELSE
    b_next_inf = b_initial2-0.075d0
    b_next_sup = b_initial2+0.075d0
ENDIF

IF (psi1.LT.500d0) THEN
    b_next_inf = bmin    
    b_next_sup = bmax
ENDIF

old_value = 10d+33
do i=1,num
   b_next_grid(i) = b_next_inf + (b_next_sup - b_next_inf) * (i-1)/ (num - 1)
   new_value = objective_fun_gN(b_next_grid(i),a_initial2,b_initial2,i_default_global2)

   if (new_value <= old_value) then
      index_opt = i
      b_next_guess = b_next_grid(i)
      old_value = new_value
  end if
   if ((DABS(old_value-new_value) <= 1d-6)) then
      index_opt2 = i
      old_value2 = new_value
   end if
end do

STEP = (b_next_grid(2) - b_next_grid(1))*0.5
BOUND = 10*EXP(amax)
XACC = 1d-6
MAXFN = 1000
   
IF (old_value.GE.1d5) THEN
    b_next = 0d0
    v_value = -old_value
ELSE        
    call DUVMIF(objective_fun2, b_next_guess, STEP, BOUND, XACC, MAXFN, b_next)
    v_value = -objective_fun2(b_next)
end if

CALL optimize_wbar(b_next2,v_value2,a_initial2,b_initial2,i_default_global2)
value3=v_value2

IF (v_value2.GT.v_value) THEN
    v_value=v_value2
    b_next=b_next2
    
    reg_global = 2
    x_global     = x_global3
    gn_global    = gn_global3
    ctrad_global = ct_global3  
ENDIF

CONTAINS

REAL(8) FUNCTION objective_fun2(x)

REAL(8),INTENT(IN)::x

objective_fun2 = objective_fun_gN(x,a_initial2,b_initial2,i_default_global2)

END FUNCTION
end subroutine
!**********************************************************************************
subroutine optimize_wbar(b_next,v_value,a_initial2,b_initial2,i_default_global2)

USE GLOBAL

implicit none

REAL(8),INTENT(IN)::a_initial2,b_initial2
INTEGER,INTENT(IN)::i_default_global2
integer :: MAXFN, i,t, d, index_opt,index_opt2, index_vector(1), num, num_tirar
parameter (num=300, num_tirar = 500)  
DOUBLE PRECISION :: b_next, v_value, objective_fun_gN_wbar, new_value, old_value,old_value2, b, vector(nb),&
                    b_next_grid(num), STEP, BOUND, XACC, b_next_guess
real(8)::b_next_inf,b_next_sup 
EXTERNAL objective_fun_gN_wbar, DUVMIF

IF (taxrate.LT.1d0) THEN
    b_next_inf = b_initial2-0.2d0
    b_next_sup = b_initial2+0.2d0
ELSE
    b_next_inf = b_initial2-0.03d0
    b_next_sup = b_initial2+0.03d0
ENDIF

IF (psi1.LT.500d0) THEN
    b_next_inf = bmin   
    b_next_sup = bmax
ENDIF

old_value = 10d+33
do i=1,num
   b_next_grid(i) = b_next_inf + (b_next_sup - b_next_inf) * (i-1)/ (num - 1)
   new_value = objective_fun_gN_wbar(b_next_grid(i),a_initial2,b_initial2,i_default_global2)
   if (new_value <= old_value) then
      index_opt = i
      b_next_guess = b_next_grid(i)
      old_value = new_value
  end if
   if (DABS(old_value-new_value) <= 1d-6) then
      index_opt2 = i
      old_value2 = new_value
  end if
end do

STEP = (b_next_grid(2) - b_next_grid(1))*0.5
BOUND = 10*EXP(amax)
XACC = 1d-6
MAXFN = 1000
   
IF (old_value.GE.1d5) THEN
    b_next = 0d0
    v_value = -old_value
ELSE        
    call DUVMIF(objective_fun2, b_next_guess, STEP, BOUND, XACC, MAXFN, b_next)
    v_value = -objective_fun2(b_next)
end if


CONTAINS

REAL(8) FUNCTION objective_fun2(x)

REAL(8),INTENT(IN)::x

objective_fun2 = objective_fun_gN_wbar(x,a_initial2,b_initial2,i_default_global2)

END FUNCTION
end subroutine
!**********************************************************************************
subroutine iterate_gN

USE GLOBAL
USE OMP_LIB

IMPLICIT NONE

INTEGER :: d

REAL(8)::a_initial,b_initial
INTEGER::i_default_global
REAL(8)::a_valor,b_valor,b_valor2,v_valor,convergence,criteria,deviation,q_fun,b_next,bmin_old
REAL(8)::b0_next,b_next1(na),b,w_value,dev_q_paid,dev_q_scheme,dev_v
REAL(8)::v0_value,v1_value,q,interpolate,b_next0,b_grid_local(nb),bmin_local
REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb),DCSVAL,epsilongrid(ncdf), dev_q_nodef

DOUBLE PRECISION, DIMENSION(nb,na)::v0_matrix_new, vaut_matrix_new, q_matrix_new, default_decision_new,q_scheme, q_scheme_new
DOUBLE PRECISION, DIMENSION(nb,na)::v_matrix_new,dev_matrix,q_matrix_nodef_new
REAL(8)::yTmax,value1
INTEGER i_b, i_a, i, j, i_b_zero, indices(1:2), indice_b, counter, ILEFT, IRIGHT,NINTV,DNORIN,new_index_b,index_b 
REAL(8)::v0_fun,x0,x1(na),b_fun,bmin_feas,objective_u_fe_gN,objective_u_unemp_gN,objective_fun_gN
REAL(8)::objective_u_fe,objective_u_unemp,objective_fun
REAL(8)::min_m,natlevel,value2,gntemp(na),cttemp(na),ct_matrix_new(nb,na)
EXTERNAL::objective_u_fe,objective_u_unemp,objective_fun
EXTERNAL::objective_u_fe_gN,objective_u_unemp_gN,objective_fun_gN
REAL(8)::new_value,x_globalx,gn_globalx,ctrad_globalx ,v_valor1,v_valor2,b_next01,b_next02
INTEGER::index_bfeas(na),index_opt,ii_b,nbaux
CHARACTER(len=29)::filename

criteria = tol   !CRITERIA FOR CONVERGENCE
convergence = -1
counter  = 0
bmin_old = bmin
vaut_matrix_new = vaut_matrix
v0_matrix_new   = v0_matrix
v_matrix_new    = v_matrix
q_scheme = 0d+0

do WHILE((convergence<0).AND.(counter.LE.maxiter))
    index_bfeas4 = 0d0
    
    !to compute spline for feasible debt gridpoints
    do i_a = 1,na
        new_value=-1d5
        i_b=0
        index_opt = 0
        
        do while (new_value <= -1d3)
           index_opt = i_b
           IF (i_b.EQ.nb) EXIT
           i_b=i_b+1
           new_value = v0_matrix(i_b,i_a)
        end do
        index_bfeas4(i_a)=index_opt     
    ENDDO

    index_bfeas2=MAXVAL(index_bfeas4)        

      deviation = 0d+0
      counter = counter + 1
      dev_v   = 0d+0
      dev_q_scheme = 0d+0
      dev_q_paid   = 0d+0
      dev_q_nodef  = 0d+0    
         
     do i_a = 1,na         
         IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
             ii_b=index_bfeas4(i_a)+1
             nbaux =nb-ii_b+1
             FDATA(1:nbaux) = v0_matrix(ii_b:nb,i_a)
	         ILEFT  = 0
	         IRIGHT = 0
	         DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(ii_b+1) - bgrid(ii_b))
	         DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	         call  DCSDEC (nbaux, bgrid(ii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, ii_b:nb), coeff_matrix(i_a, :,ii_b:nb))
         ENDIF         
     ENDDO
     
        do i_a = 1,na
            FDATA = q_matrix_nodef(:,i_a)
	        ILEFT  = 0
	        IRIGHT = 0
            DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(2) - bgrid(1))
	        DRIGHT = (FDATA(nb)-FDATA(nb-1)) / (bgrid(nb) - bgrid(nb-1))
	        call  DCSDEC (nb, bgrid, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, BREAK_GRID, CSCOEF)
            break_matrix_q(i_a, :) = BREAK_GRID
            coeff_matrix_q(i_a, :,:) = CSCOEF
        end do
        
    IF (psi1.LE.500d0) THEN
    !! Compute the value of repayment
    !!!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(a_initial,b_initial,i_default_global,b_valor,v_valor)
    !!!$OMP DO COLLAPSE(1)
    LOOP_A1: do i_a = 1,na
              
                IF (choice_constant_gn.EQ.1)  constant_gn = gnaut_base(i_a)
              
                 a_initial = agrid(i_a)
                 b_initial = bgrid(1)
    
                  i_default_global=2 !Country defaults
                  b_valor = 0d0
                  IF (only_repay.EQ.1) THEN
                      v_valor= -objective_fun(b_valor,a_initial,b_initial,i_default_global)
                  ELSEIF (only_repay.EQ.0) THEN
                      v_valor= -objective_fun_gN(b_valor,a_initial,b_initial,i_default_global)
                  ENDIF                  
                  x1(i_a)     = x_global
                  gntemp(i_a) = gn_global
                  cttemp(i_a) = ctrad_global
                  b_next1(i_a) = b_valor
                  vaut_matrix_new(:, i_a) = v_valor
            end do LOOP_A1
            !!!$OMP END DO
            !!!$OMP END PARALLEL
            ELSE
                  x1      = 0d0
                  gntemp  = 0d0
                  cttemp  = 0d0
                  b_next1 = 0d0
                  vaut_matrix_new = -1d10
            ENDIF
    
    !!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(a_initial,b_initial,i_default_global,v_valor,b_next0)
    !!$OMP DO COLLAPSE(2)
    LOOP_A2: do i_a = 1,na
    LOOP_B:  do i_b = 1,nb 
                  a_initial = agrid(i_a)
                  b_initial = bgrid(i_b)

                  IF (choice_constant_gn.EQ.1)  constant_gn = constant_gn_value+coefyt*a_initial
                  i_default_global=1 !Country does not default
                  
                  call optimize_gN(b_next02, v_valor2,a_initial,b_initial,i_default_global)
                  
                   
                  b_next0 = b_next02
                  adhere_rule(i_b, i_a)=1
                  v_valor = v_valor2

                  v0_matrix_new(i_b, i_a) = v_valor
                  q_matrix_nodef_new(i_b, i_a) = q_fun(b_next0, a_initial)
                  
                  x0 = x_global

                  if (vaut_matrix_new(i_b, i_a) > v0_matrix_new(i_b, i_a)) THEN 
                     b_next_matrix(i_b, i_a) = 0d0 
                     x_matrix_new(i_b, i_a)  = x1(i_a)                     
                     gn_matrix_new(i_b, i_a) = gntemp(i_a)     
                     ct_matrix_new(i_b, i_a) = cttemp(i_a)   
                     default_decision_new(i_b,i_a) = 2
                     v_matrix_new(i_b, i_a) = vaut_matrix_new(i_b, i_a)
                  else
                    b_next_matrix(i_b, i_a) = b_next0                     
                    gn_matrix_new(i_b, i_a) = gn_global   
                    x_matrix_new(i_b, i_a)  = x0                     
                    ct_matrix_new(i_b, i_a) = ctrad_global   
                        
                    default_decision_new(i_b, i_a) = 1
                    v_matrix_new(i_b, i_a) = v0_matrix_new(i_b, i_a)
                  END if
                  
                        q_scheme_new(i_b, i_a) =  q_fun(b_initial,a_initial)
                        q_matrix_new(i_b, i_a) = q_fun(b_next_matrix(i_b,i_a),a_initial)
                  
                   dev_v =        max(ABS(v_matrix_new(i_b, i_a) - v_matrix(i_b, i_a)), dev_v)
                   dev_q_paid =   MAX(ABS(q_matrix_new(i_b, i_a) - q_matrix(i_b, i_a)), dev_q_paid)
                   dev_q_scheme = MAX(ABS(q_scheme_new(i_b, i_a) - q_scheme(i_b, i_a)), dev_q_scheme)
                   dev_q_nodef =   MAX(ABS(q_matrix_nodef_new(i_b, i_a) - q_matrix_nodef(i_b, i_a)), dev_q_nodef)

                   deviation = MAX(deviation, MAX(dev_v, dev_q_scheme))
                   dev_matrix(i_b, i_a) = (q_matrix_new( i_b, i_a) - q_matrix( i_b, i_a))

    end do LOOP_B
     end do LOOP_A2
    !!!$OMP END DO
    !!!$OMP END PARALLEL           
                

 WRITE(*, '(I5, X, F15.6, X, F15.6, X, F12.3, X, F12.3)') counter, dev_v,dev_q_scheme,MINVAL(-b_next_matrix*(v0_matrix_new.GE.-500d0)),MAXVAL(b_next_matrix)

    
IF ((save_results.EQ.1)) THEN
      open (10, FILE='graphs/v.txt',STATUS='replace')
      open (11, FILE='graphs/default.txt',STATUS='replace')
      open (12, FILE='graphs/q.txt', STATUS='replace')
      open (13, FILE='graphs/b_next.txt',STATUS='replace')
      open (14, FILE='graphs/devq.txt',STATUS='replace')
      open (15, FILE='graphs/dev_iter.txt',ACCESS='append')
      open (16, FILE='graphs/q_paid.txt', STATUS='replace')
      open (110, FILE='graphs/b_grid.txt',STATUS='replace')
      open (113, FILE='graphs/bounds.txt',STATUS='replace')
      open (114, FILE='graphs/x.txt',STATUS='replace')
      open (115, FILE='graphs/gn.txt',STATUS='replace')
      open (116, FILE='graphs/ctt.txt',STATUS='replace')
      
      WRITE(113, '(F15.11, X, F15.11)') bmin, bmax
      WRITE(113, '(F15.11, X, F15.11)') amin, amax

      WRITE(15,'(I3,X,F15.10, X, F15.10)') counter,dev_v ,dev_q_paid

         do i_b = 1,nb
            WRITE(110, '(F15.11)') bgrid(i_b)
            do i_a = 1,na
                     a_initial = agrid(i_a)
                     b_initial = bgrid(i_b)

                     WRITE(10, '(F18.10, X, F18.10, X, F18.10)') MAX(-1000d0,v_matrix(i_b, i_a)), &
                     MAX(-1000d0,v0_matrix(i_b, i_a)), MAX(-1000d0,vaut_matrix(i_b, i_a))
                     WRITE(11, '(F6.2)') defaultgrid(default_decision(i_b, i_a))

                     WRITE(12, '(F15.11)') q_fun(bgrid(i_b), agrid(i_a))

                     WRITE(13, '(F15.11)') b_next_matrix(i_b, i_a)
                     WRITE(14, '(F15.11)') dev_matrix(i_b, i_a)
                     WRITE(16, '(F15.11, X, F15.11)') q_matrix_new(i_b, i_a), q_matrix_nodef_new(i_b, i_a)
                     WRITE(114, '(F15.11)') x_matrix_new(i_b, i_a)
                     WRITE(115, '(F15.11)') gn_matrix_new(i_b, i_a)                     
                     WRITE(116, '(F15.11)') ct_matrix_new(i_b, i_a)                     
            end do
         end do

 CLOSE(10)
 CLOSE(11)
 CLOSE(12)
 CLOSE(13)
 CLOSE(14)
 CLOSE(15)
 CLOSE(16)
 CLOSE(110)
 CLOSE(113)
 CLOSE(114)
 CLOSE(115)
 CLOSE(116)

ENDIF

!UPDATE VALUES OF MATRICES
   default_decision = default_decision_new
   if (deviation < criteria) then
      convergence =1
   end if
    

       q_matrix = q_matrix_new
       q_scheme = q_scheme_new
       q_matrix_nodef = q_matrix_nodef_new
       vaut_matrix = vaut_matrix_new
       v0_matrix = v0_matrix_new
       v_matrix = max(v0_matrix, vaut_matrix)
    IF (counter.EQ.maxiter) EXIT

end do

counter_global = counter
 
PRINT*, 'Model solved'

end subroutine iterate_gN
!************************************************************************************************************
SUBROUTINE pol_fcns_gN
USE GLOBAL

IMPLICIT NONE

REAL(8)::a_initial,b_initial,a0,aux0
INTEGER::i_default_global
REAL(8)::a,b
INTEGER::i_a,i_b,i_bp,i,j,i0(1),zero_h,zero_h2
INTEGER,PARAMETER::nb3=501

REAL(8)::wmax,p0max,h0max,gn0max,ct0max,tau0max,ww0max,q_fun
REAL(8)::objective_fun_gN,objective_fun,interpolate,v_valor,b_next
REAL(8)::ctval,cnval,p0,q,ca1(na2,nb2),v_valor1,v_valor2,b_next01,b_next02
EXTERNAL::q_fun,objective_fun_gN,objective_fun,interpolate
REAL(8)::bgrid_base(nb3),h_baseX(na,nb3),gn_baseX(na,nb3)    
REAL(8)::gn_base2(na2,nb2),h_base2(na2,nb2),gn_globalx,x_globalx,ctrad_globalx   
 
    
do i=1,na2
   agrid2(i) = amin + (amax - amin) * REAL((i-1d0) / (na2 - 1d0))
end do

do i=1,nb2
   bgrid2(i) = bmin + (bmax - bmin) * REAL((i-1d0) / (nb2 - 1d0))
end do

i0 = LOCATE2(bgrid2,0d0)
i_bzero2 = i0(1)
bgrid2(i_bzero2)=0d0
choice_glob=1

OPEN(1,FILE='output//bgrid2.txt')           
DO i = 1,nb2
    WRITE(1,*) bgrid2(i)
ENDDO
CLOSE(1)

OPEN(1,FILE='output//agrid2.txt')           
DO i = 1,na2
    WRITE(1,*) agrid2(i)
ENDDO    
CLOSE(1)
      
CALL BUILD_TRANSITION(na2,agrid2,rhoa,sigmaa,mua,prob)
   
DO i_a=1,na2

    a_initial = agrid2(i_a)
    a0 = DEXP(a_initial)

    b_initial = 0d0 
    i_default_global = 2  !default occurs
    b = 0d0

    IF (choice_constant_gn.EQ.1)  constant_gn = constant_gn_value
                    
    IF (only_repay.EQ.1) THEN
        v_valor = -objective_fun(b,a_initial,b_initial,i_default_global)
    ELSEIF (only_repay.EQ.0) THEN
        v_valor = -objective_fun_gN(b,a_initial,b_initial,i_default_global)
    ENDIF
    
    vaut1(i_a)  = v_valor
    ctval  = ctrad_global 
    cnval  = x_global**alpha-gn_global
    p0     = (1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)
    
    paut1(i_a)  = p0
    haut1(i_a)  = x_global
    gnaut1(i_a) = gn_global
    ctaut1(i_a) = ctval
    tauaut1(i_a)= taxrate*(a0+in*p0*x_global**alpha)
    waut1(i_a)  = alpha*p0*x_global**(alpha-1d0)
    
    DO i_b=1,nb2    
            
    i_default_global = 1
    b_initial = bgrid2(i_b)
    
    b = interpolate(b_initial, a_initial, b_next_matrix)
    i_bp = LOCATE2(bgrid2,b)
    
    IF (choice_constant_gn.EQ.1)  constant_gn = constant_gn_value+coefyt*a_initial

    call optimize_gN(b_next02, v_valor2,a_initial,b_initial,i_default_global)
                  
    b_next = b_next02
    adhere_rule1(i_a,i_b)=1
    v_valor = v_valor2
    vcon1(i_a,i_b)  = v_valor

    
    bnext1(i_a,i_b) = b_next
    q = q_fun(b_next,a_initial)

    aux0 = -q*(b_next-(1d0-lambda)*b_initial ) + payoff*b_initial            !p*gn - T

    ctval  = ctrad_global
    cnval  = x_global**alpha-gn_global
    p0     = (1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)    
    bp1(i_a,i_b)    = LOCATE2(bgrid2,b_next)
    p1(i_a,i_b)     = p0
    h1(i_a,i_b)     = x_global
    gn1(i_a,i_b)    = gn_global
    ct1(i_a,i_b)    = ctval
    q1(i_a,i_b) = q_fun(b_initial,a_initial)    

    ca1(i_a,i_b)    = aux0
    tau1(i_a,i_b)   = taxrate*(a0+in*p0*x_global**alpha)
    w1(i_a,i_b)     = alpha*p0*x_global**(alpha-1d0)   

    def1(i_a,i_b) = 1d0
    IF (vaut1(i_a).GT.vcon1(i_a,i_b)) def1(i_a,i_b) = 0d0
ENDDO
ENDDO

OPEN(1,FILE='output//gvcon.txt')
OPEN(2,FILE='output//gp.txt')
OPEN(3,FILE='output//gh.txt')
OPEN(4,FILE='output//gibp.txt')
OPEN(5,FILE='output//gct.txt')
OPEN(6,FILE='output//gbnext.txt')
OPEN(7,FILE='output//ggn.txt')
OPEN(8,FILE='output//gtau.txt')    
OPEN(9,FILE='output//gw.txt')    
OPEN(10,FILE='output//gq.txt')    
OPEN(11,FILE='output//gdef1.txt')    
OPEN(12,FILE='output//gca.txt')    
OPEN(13,FILE='output//gadhere_rule.txt')    
            
DO i = 1,na2
DO j = 1,nb2
    WRITE(1,*) vcon1(i,j)
    WRITE(2,*) p1(i,j)
    WRITE(3,*) h1(i,j)
    WRITE(4,'(I7)') bp1(i,j)
    WRITE(6,*) bnext1(i,j)
    WRITE(5,*) ct1(i,j)
    WRITE(7,*) gn1(i,j)
    WRITE(8,*) tau1(i,j)
    WRITE(9,*) w1(i,j)
    WRITE(10,*) q1(i,j)
    WRITE(11,*) def1(i,j)
    WRITE(12,*) ca1(i,j)
    WRITE(13,'(I7)') adhere_rule1(i,j)
ENDDO
ENDDO
    
CLOSE(1)
CLOSE(2)
CLOSE(3)
CLOSE(4)
CLOSE(5)
CLOSE(6)
CLOSE(7)
CLOSE(8)
CLOSE(9)
CLOSE(10)
CLOSE(11)
CLOSE(12)
CLOSE(13)
    
OPEN(1,FILE='output//gvaut.txt')
OPEN(2,FILE='output//gpaut.txt')
OPEN(3,FILE='output//ghaut.txt')
OPEN(5,FILE='output//gctaut.txt')
OPEN(7,FILE='output//ggnaut.txt')
OPEN(8,FILE='output//gtauaut.txt')    
OPEN(9,FILE='output//gwaut.txt')    
            
DO i = 1,na2
    WRITE(1,*) vaut1(i)
    WRITE(2,*) paut1(i)
    WRITE(3,*) haut1(i)
    WRITE(5,*) ctaut1(i)
    WRITE(7,*) gnaut1(i)
    WRITE(8,*) tauaut1(i)
    WRITE(9,*) waut1(i)
ENDDO
    
CLOSE(1)
CLOSE(2)
CLOSE(3)
CLOSE(5)
CLOSE(7)
CLOSE(8)
CLOSE(9)
    
END SUBROUTINE pol_fcns_gN
!****************************************************************************************